Zitterbewegung of Klein-Gordon particles and its simulation by classical systems 
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The Klein-Gordon equation is used to calculate the Zitterbewegung (ZB, trembling motion) of 
spin-zero particles in absence of fields and in the presence of an external magnetic field. Both 
Hamiltonian and wave formalisms are employed to describe ZB and their results are compared. It 
is demonstrated that, if one uses wave packets to represent particles, the ZB motion has a decaying 
behavior. It is also shown that the trembling motion is caused by an interference of two sub- 
packets composed of positive and negative energy states which propagate with different velocities. 
In the presence of a magnetic field the quantization of energy spectrum results in many interband 
frequencies contributing to ZB oscillations and the motion follows a collapse-revival pattern. In the 
limit of non-relativistic velocities the interband ZB components vanish and the motion is reduced to 
cyclotron oscillations. The exact dynamics of a charged Klein-Gordon particle in the presence of a 
magnetic field is described on an operator level. The trembling motion of a KG particle in absence 
of fields is simulated using a classical model proposed by Morse and Feshbach - it is shown that a 
variance of a Gaussian wave packet exhibits ZB oscillations. 

PACS numbers: 03.65.Pm, 11.40.-q, 03.65.-w 



INTRODUCTION 



The phenomenon of Zitterbewegung (ZB, trembling 
motion) goes back to Schrodinger who proposed it in 1930 
for free relativistic electrons in a vacuum [l| . Schrodinger 
observed that, due to non-commutativity of the velocity 
operators with the Dirac Hamiltonian, relativistic elec- 
trons experience a trembling motion in absence of ex- 
ternal fields. It was later recognized that ZB is due to 
an interference of electron states with positive and neg- 
ative electron energies. A very high frequency of ZB 
in a vacuum, corresponding to Hujz = 2m e c 2 , and its 
very small amplitude on the order of the Compton wave- 
length A c = h/m e c ~ 3.86 x 10 -3 A made it impossible to 
observe this effect in its original form with the currently 
available experimental methods. However, in a recent 
work Gerritsma et al. Q simulated the 1+1 Dirac equa- 
tion and the resulting Zitterbewegung with the use of 
trapped ions excited by laser beams. The important ad- 
vantage of this method is that one can simulate also the 
basic parameters of the Dirac equation and tailor their 
desired values. The result of Gerritsma et al. allows one 
to expect that observable effects for relativistic particles 
in a vacuum can be convincingly reproduced with more 
"user friendly" parameters. In general, there has been 
recently a revival of interest in the relativistic-type equa- 
tions related to "the rise of graphene" Q, topological 
insulators and similar systems in narrow-gap semicon- 
ductors Q. 

The purpose of our paper is to describe the phe- 
nomenon of Zitterbewegung for charged Klein- Gordon 
(KG) spin-zero particles in absence of fields and in the 
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presence of a magnetic field 0-0] • The Zitterbewegung 
of KG particles in absence of fields was described be- 
fore, see |3-fl0l|. However, in our treatment we introduce 
a number of additional elements. First, we describe the 
particles by wave packets and show that this feature leads 
to a transient character of the resulting ZB motion. Sec- 
ond, we use both the Hamiltonian and wave forms of 
Klein-Gordon equation (KGE) and show the equivalence 
of the two approaches. Third, we point out that ZB 
is a result of interference between positive and negative 
energy sub-packets propagating with different velocities. 
Fourth, we simulate classically the ZB motion using a 
simple mechanical system proposed by Morse and Fesh- 
bach 11]. Still, our main objective is to consider in detail 
the dynamics of a charged KG particle in the presence 
of an external uniform magnetic field and describe the 
phenomenon of ZB in this situation. To the best of our 
knowledge this problem has not been treated before. 

The one-particle Klein-Gordon equation for spin-zero 
particles leads to some well known difficulties fiol . [l^ . 
The KG equation involves second time derivative, the 
probability density is not positively definite, there are 
problems with the position operator or vanishing square 
of the velocity operator. For this reason in the present 
work we calculate ZB of average current which has well 
defined meaning in the theory of KG equation. For 
charged particles the average current is proportional to 
average particle velocity, so in our work we calculate one 
of these two quantities. In previous treatments of ZB 
for Dirac equation, simulation by trapped ions or solid- 
state systems, the authors usually calculated ZB of the 
position operator. 

In our considerations we encounter another interesting 
anomaly of KG equation, namely, that particle velocities 
can exceed the speed of light for sufficiently large mo- 
menta. In other words it appears that, in contrast to 
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the Dirac equation for electrons, KGE does not posses 
an automatic "safety brake" for velocities to keep them 
below c. To our knowledge this feature has not been 
remarked before, so we mention it throughout our work. 

Our paper is organized as follows. In Section[H]we cal- 
culate ZB of a wave packet using the Hamiltonian formal- 
ism, in Section IlIII we obtain similar results with the use 
of KG waves and discuss explicitly physical background 
for the transient behavior of ZB motion. Section |IV] con- 
tains a description of ZB for a charged KG particle in a 
magnetic field. In Section |V] we simulate classically the 
ZB phenomenon using a system proposed by Morse and 
Feshbach. In Section IVII we discuss our results, the pa- 
per is concluded by a summary. Appendix [A] contains a 
derivation of particle dynamics in the presence of a mag- 
netic field, Appendix [B] discusses the problem of high 
particle velocities, in Appendices ICl and IDl we give some 
mathematical details. 



II. ZITTERBEWEGUNG IN VACUUM 

We begin by considering a Klein-Gordon particle in 
absence of external fields. The Klein-Gordon equation in 
the Hamiltonian form is [lJ 



<9* 

ih— = H^>. 

dt 



Here the Hamiltonian is 



T 3 + IT-2 „ 2 2 

H = — p + r 3 mc , 

Am 



(1) 



(2) 



where m is particle mass, p is particle momentum 
and Tj (j = 1,2,3) are the Pauli matrices Oj, respec- 
tively. The wave function ^ is a two-component vector 



* = 



(3) 



In the Hamiltonian form one can introduce the Heisen- 
berg picture f]~3j j . The z-th component of the time- 
dependent velocity operator is 



v z (t) = e lIit ' h v z {Q)e~ l " t/r \ 



(4) 



where v z (0) = dH/dp z . In this representation v z (t) is 
a 2 x 2 matrix operator. Expanding e lHt / h = 1 + Ht + 
(l/2!)# 2 + . . . and noting that H 2 = E 2 , where the en- 
ergy is E = ±cpo with 



Po 



= +y m 2 c 2 + p 2 , 



we obtain 



jHt/h 



iH 

cos(Et/H) + — sm(Et/H) 
E 



(5) 



(6) 



The velocity operator in Eq. ^ is a product of three 
matrices. Its (1, 1) component is 

(«,) u (t) = h + |% [cos(2Et/h) - 1] . (7) 



The remaining elements of v z (t) are calculated similarly. 
The v x and v y components of the velocity operator are 
obtained from v z (t) by the replacement p z — > p X7 p y , re- 
spectively. In the non-relativistic limit p -C mc we obtain 
in Eq. (J7J) the classical motion (v z )n(t) ~ p z /m. In ab- 
sence of external fields p% are good quantum numbers. 
We introduce p = hk and q — \ c k, where the effective 
Compton wavelength is A c = h/mc. Also, we introduce 
a useful frequency ujq = (mc 2 )/H. Both A c and wo refer 
to particles of mass m. In the above notation Eq. ([7]) 
becomes 



(v z )ll(t) = cq z + 



2(i + g 2 ) 



3(2w tVl + q 2 ) - 1 



(8) 

The first term in Eq. (O corresponds to the classi- 
cal motion of a particle while the second term describes 
rapid oscillations of the velocity. The velocity oscillates 
from v max = cq z to v min — cq z /(l + q 2 ). Since the max- 
imum velocity of the particle is c, there must be \q\ < 1. 
We notice that, in principle, Eq. ((5J admits velocities 
above the speed of light. We discuss this issue in more 
detail in Appendix [B] The frequency of oscillations varies 
from lu — 2luq for low q to cj — 2\[2ujq for \q\ — 1. The 
velocity oscillations taking place in absence of external 
fields are called Zitterbewegung. 

Integrating (v z )n(t) in Eq. (jSJ over time we have 



zii(t) 



/r\\ , , c q J -q z 
zn(0) + cq z t - -— — -t 



2 1 



q 2 q z 



4 (l + 9 2 )3/2 



Sill 



(2w tv/r 



(9) 



The amplitude of ZB oscillations of the position operator 
is on the order of A c . The operator Zn(t) is obtained 
in the formal way, physical limitations to the position 
operator will be discussed below. 

In order to obtain physical observables one needs to av- 
erage the operator quantities over the wave packet. The 
average velocity (v z (t)) of the wave packet \W) is 

(v z (t)) = (W\r 3 v z (t)\W) =J2( W \p) ( P\ t ^Mp')¥\W). 

pp' 

( 10 ) 

For KGE in the Hamiltonian form the matrix elements 
of operators include an additional t 3 factor [l3| . We take 
the wave packet in the form of a two-component vec- 
tor (r\W) = (l,0) T (r|w) with one non- vanishing compo- 
nent. Here (r\w) = w(r) is a Gaussian function with a 
nonzero momentum hko 



w(r) 



1 



■cxp[-r 2 /(2d 2 ) +ik Q r}. (11) 



(<V^) 3 / 2 

There is w{k) = j e~ lkr / h w(r)d 3 r and we have 

(k\w) = (2dv^F) 3/2 exp[-d 2 (fe - fc ) 2 /2]. (12) 
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duce the unity operator Q 

i = ^2 \ks)(ks\r 3 s, 

ks 

where s = ±1, and 
(r\ks) 



e / mc + spo 



2^/mcp^ \ mc- sp 



(14) 



(15) 



are the two eigenstates of H corresponding to the posi- 
tive and negative energies E s = scpo. These states are 
normalized according to (fcs|T 3 |fcV)/(27r) 3 / 2 = s5kk>5 ss >. 
Then 

\W) =Y J s\ks){ks\n\W) = J2s\ks)W ks , (16) 
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ks 
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FIG. 1: Calculated velocity of wave packet in absence of ex- 
ternal fields for three packet widths d. The phenomenon of 
transient Zitterbewegung is seen. The initial packet wave vec- 
tor is feo = (0,0, koz) with ko z = 0.8A^~ . Time is expressed 
in t c = h/mc 2 units. Initial packet velocity is vo z = hkoz/m — 
0.8c, its final velocity depends on packet parameters. 



The wave packet \W) selects (1,1) component of the ve- 
locity matrix v z (t). From Eqs. (|10|) and (fl~2j) we obtain 



where W ks = (ks\r 3 \W). The sub-packet of pos- 
itive energy states is \W+) = J2k \k+)W k +, while 
the sub-packet of negative energy states is \W— ) = 
Efe \k-)W k — Using Eqs. JT5J) and (JIBJ) we find 

TT fes = {2d^f/ ^ mC + SPa) e-^-^ 2 l\ (17) 
2^/mcp 

The average packet velocity is 

(«.(*)> = E ss'W^ s W kls/ (ks\T 3 v z (t)\k's') 



kk' ss' 



dH , 



V S s'W^W k , s ,e^-^^{k S \r 3 — \k' S '). (18) 

^ OPz 



(«.(*)> 



-^2 y exp[-d 2 (g-g ) 2 



1 f? z 



2 1 



i(u> a ty/T 



We defined w s = su)q^/1 + (kX c ) 2 and used the equality 
d 3 q, (13) (fes|r 3 e i#t / ft = (fcs|e^ tt/n T 3 = e^fcshi. (19) 



where <i c = d/X c . This integral is nonzero only if q has 
a nonzero z-th component, so we take qo — (0,0, qoz) T ■ 
Selecting the z axis to be parallel to qo and using the 
spherical coordinates we calculate the integrals over the 
two angular variables. The remaining integral over q is 
computed numerically. 

In Fig. Q] we plot the average packet velocity (v z (t)) cal- 
culated from Eq. (fi"3")) for three different packet widths d. 
The time on the horizontal axis is expressed in t c — 
H/mc 2 units, where t c = (m e /m) x 1.29 x 10~ 21 s and m e 
is the electron mass. In all cases the motion has a tran- 
sient character. Physically, the decay of ZB oscillations 
is due to different propagation velocities of sub-packets 
corresponding to the positive and negative energy states. 
We analyze this effect below. It is seen that the final 
packet velocity differs from the initial value hko z /m. In 
the limit of d — > oo the velocity oscillations do not decay 
in time. 

Now we calculate the average velocity by splitting the 
initial wave packet into two sub-packets corresponding to 
the positive and negative energy states. First we intro- 



which follows from the properties: H = t 3 Wt 3 
and (fcs|#t = (H\ks)Y = E s (ks\. Another proof of the 
identity (|T9|) is given in Appendix [C] There is also 

(fe s |r 3 |^|feV) - (2^) 3 ^<W, (20) 

which does not depend on s and s' . Combining Eqs. (|18[) 
- (l20l) we obtain 



(21) 



(27r) d m J pfi 



ss' {mc + spo)(mc + s'p^e 



The average velocity in Eq. ([21]) is a sum of four terms. 
The term with s = s' = +1 describes the motion of 
positive energy sub-packet, while the term with s = s' = 
— 1 corresponds to the negative energy sub-packet 



d 3 



4m7r 3 / 2 



J U ± Pz e- d2 ( k - k ^d 3 k. (22) 
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Thus the two sub-packets move with different velocities. 
Their relative velocity is 



\ rel 



cd 3 I Vz_ e -£ 

7T 3/2 J PO 



(fc-fc ) 2 d 3 fc> ( 23 ) 



Two terms in Eq. (|2ip with s s , corresponding to an 
interference of the two packets, give rise to an oscillatory 
term 



4m7r 3 / 2 ./I Po ^ 



x cos(2uj k t)e- d2{k - ko '> 2 d 3 k, (24) 



where uj^ — y/T+ (kX c ) 2 . According to the Riemann- 
Les beg ues theorem this term has a transient charac- 
ter pjj]. Performing integrations in Eqs. (f22j) and (|24|) . we 
obtain again Eq. f| 13[) . Thus we showed that the ZB oscil- 
lations arise from the interference of positive and negative 
energy states. After a certain time the two sub-packets 
are sufficiently far away from each other and the overlap 
between them vanishes, which results in the disappear- 
ance of ZB oscillations. This explains the behavior of 
velocity shown in Fig. [TJ 

To evaluate the decay of ZB oscillations, we estimate 
the time after which the two sub-packets will be sep- 
arated from each other by the distance 2c?. Assuming 
that fcoA c ~ 1, the relative velocity between the two 
sub-packets is (v z ) rel ~ c(ko\ c ). The time interval after 
which the distance between the sub-packets exceeds 2c? is 



2d 
ck X c 



(25) 



It is seen in Fig. [JJthat the ZB oscillations nearly disap- 
pear after td- For example there is td — 5t c for d = 2X C . 
Since the ZB frequency is 2ujq = 2mc 2 /h, a number of 
non- vanishing oscillations is approximately 



N„. 



2uotd 
2ir 



1 



k X c 



(26) 



The above estimation correctly evaluates the number of 
ZB oscillations seen in Fig. [TJ The optimal conditions for 
an appearance of ZB are: wide packets and small values 
of |feo|. On the other hand, for too small values of |feo| 
one of the two sub-packets disappears, see Eq. (|22[), which 
reduces amplitude of ZB oscillations. 



III. WAVE FORM OF KGE 

Now we intend to demonstrate a relation between the 
ZB oscillations of the average packet velocity calculated 
above with the use of the Hamiltonian form of KGE and 
an average current obtained from the wave form of KGE. 
In absence of external fields the Klein-Gordon equation 
has the wave equation form 



cj>{x) - V 2 (j)(x) + 



2 2 

m c 



<t>{x) = 0, (27) 



where x — (ct,r) is the position four- vector [Ioj |. The 
solution of this equation is 



4>(x) 



(2 



2_ [ .[^\a(k)e- ik - x + b*(k)e ik - x ]d 3 k, 

7Tj J J V PO 



(28) 



where k = (ui k /c,k), uj k = w ^l + (fcA c ) 2 , 
and a(k), b*(k) are complex coefficients. Function tj> is 
normalized to 
ih 



2mc 2 



d A r = Q, (29) 



where Q = ±1 for charged particles and Q — for neutral 
particles. In the following we select Q = +1, which leads 
to 



d 3 k [a*(k)a(k) - b*(k)b(k)} 



1. 



(30) 



To determine the coefficients a(k) and b*(k) we need two 
boundary conditions for <fi and d<fi/dt at x = (0, r). Hav- 
ing specified a(k) and b*(k) one can calculate the current 
density j(x) 

j(^) = ^-[r(v0)-(v^], (si) 

and the average current (j(t)) = J j(x)d 3 r. 

Our aim is to find a correspondence between the aver- 
age packet velocity calculated in Eq. (fT"3|) and the average 
current (j(t)} given in Eq. (|3~T1) . To this end we select the 
coefficients a(k) and b*(k) in such a way that the func- 
tion (f> in the wave form of KGE corresponds to the wave 
packet (w(r),0) T in the Hamiltonian form of KGE. Re- 
lations between </>, dip/dt and the two-component wave 
function = (ip, x) T in the Hamiltonian form of KGE 
are 

4> = f + X, (32) 
id(j)/dt = mc 2 (tp-x)lft- (33) 

Since (ip, x) T = (w(r),0) T we find the coefficients a(k) 
and b*(k) from Eqs. ([5 ^1 - ([3"3" | by setting <p(t = 0,r) = 
w(r) and \ = 0. From Eq. (|3"!2)) we have 



f % f^E Uk)e +lk r + b*(k)e- lk - r ] d 3 k 
J V Po 

= (2dV^) 3/2 J e -( fc - fc «) 2d2 /2+^- d 3 fcj 
while from Eq. (|33| we have 



(34) 




-\-ik-r 



-b*(k)e- ikr ] Po d 3 k 



(2dv^) 3/2 mc / e -(k-k r^/2+ik-r d 3 k _ (35) 



d 3 ke -d 2 (k-k y/2+ik-r 



In terms including b*(k) we replace k — > — k, solve equa- 
tions (|3"4")l and ([33]) for a(k) and b*(— k), and obtain 

(2d0F)3/2 

<p(r,t) 



2(2tt) 3 

mc 
1 + 

Po 



— iufct 



1 

PO 



(36) 
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FIG. 2: Time evolution of the wave packet (absolute value) 
according to one-dimensional version of Eq. (|36[) . Initial wave 
packet (thick line) splits into two sub-packets moving with 
different velocities. Thin lines show shapes of sub-packets in 
successive time intervals 2t c = 2fi/(mc 2 ). 



The above function 4> includes both positive and negative 
energy amplitudes. For p — > there is 1 + mc/po ~ 2 
and 1 — (mc/po) ~ p 2 /2(mc) 2 . Thus the second term in 
Eq. (|36l) is much smaller than the first. In this limit the 
packet consists of the positive energy states alone. 

In Fig. [5] we plot the time evolution of the wave 
packet 4> in one dimension. The packet propagates ac- 
cording to a one-dimensional version of Eq. (|36|) . The 
initial packet is assumed in a Gaussian form 



(x,0) = —=e 

\/7T 



-x 2 /(2d 2 )+ik a x 



(37) 



Its absolute value is indicated in Fig. [2] by the thick line. 
Each thin line describes \<f>(x, t)\ in successive time inter- 
vals 2t c — 2h/ (mc 2 ). It is seen that the packet splits into 
two sub-packets moving with different velocities. The 
sub-packet at the right corresponds to positive energies 
while the sub-packet at the left corresponds to negative 
energies. The difference in the amplitudes of sub-packets 
results from different contributions of the positive and 
negative energy states in the initial packet at t = 0, see 
Eqs. (|34|) - (|35|) . The Zitterbewegung occurs only when 
the sub-packets overlap. Each of the sub-packets slowly 
spreads in time, but the spreading time is much larger 
than the overlapping time, so the ZB vanishes much 
faster than the spreading of sub-packets. 

Now we continue the calculation of average current 
given in Eq. (f3Tj) using function </> of Eq. (f36)) . This func- 
tion has the form of an integral over k. To calculate 
the spatial derivative V</> we change the order of inte- 
gration and differentiation, which can be done for any 



function decaying exponentially for k — > oo. Using the 
identity: 1 + (mc/p^) 2 = 2 — (p/po) 2 , we obtain for the 
first term of the average current 



h 

2im 



dz 



d A r 



d 3 H 



8im7r 3 / 2 



(ik z ) x 



Po 



-iuikt 



1 I e 

Po 



= I ^ lh -""{^ + 1£ WW) - «} -(38) 

Calculation of the second term in the cur- 
rent: h/(2im) f(d<fi* /dz)4>(fir : gives the same result 
but with an opposite sign, so that both terms in Eq. (|3ip 
add together. Comparing Eq. (|38]l with Eqs. (|3Tj) 
and (fll?l) we conclude that the current density (j z (t)} 
averaged over the packet (f)(x) in Eq. (|2"51) equals to the 
average velocity (v z (t)} of the packet in the Hamiltonian 
form of KGE multiplied by the particle charge. This 
way we establish an equivalence of Zitterbewegung in 
the Hamiltonian and wave equation formalisms. 

The above equivalence is valid for the average val- 
ues only. In the Hamiltonian form of KGE one can 
define the time dependent velocity operator v(t) — 
e lHt / h v(0)e~ lHt / h , which can be expressed in a closed 
form without specifying of the wave packet, see Eq. (|7|)- 
But an analogous current operator in the wave form of 
KGE can be defined as a current density j(x), which 
strongly depends on the form of function <f>. 

Even more significant differences between the Hamilto- 
nian and wave descriptions of ZB appear in the analysis 
of the position operator r(t). In the Hamiltonian form 
of KGE the position operator written in the Heisenberg 
picture is f(t) = e lHt / h f(Q) e ~ lHt / h anc f ) f or the field-free 
KGE, it can be calculated in a compact form, see Eq. © 
and Ref. ||- On the other hand, there is no well defined 
position operator r for the wave form of KGE since this 
operator is not hermitian, see Ref. [l2j]. However, one can 
calculate an average position operator for the wave form 
of KGE by integrating the average current over time 



(r(t)) = (r(0)) 



1 

Q 



[j(t))dt, 



(39) 



where the charge Q 7^ 0. This example indicates that the 
equivalence between the Zitterbewegung for the Hamilto- 
nian and wave equation formalisms holds for the average 
values only. 



IV. ZITTERBEWEGUNG IN A MAGNETIC 
FIELD 

In the presence of a magnetic field the KG Hamiltonian 
for a charged particle reads [To| 



H 



t 3 + ir 2 
2m 



(p - qA) 2 + T 3 mc 2 



(40) 
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where q is the particle charge and A is the vector poten- 
tial of a magnetic field. We assume the magnetic field B 
to be parallel to the z axis and describe it by the asym- 
metric gauge A = B(—y, 0, 0). Eigenstates of the Hamil- 
tonian arc of the form 



*(r) = e 



ik x x+ik z z 



(41) 



and the resulting eigenenergy equation is — Ety with 
H = (r 3 +iT 2 )^ [(hk x + qByf + h 2 k 2 y + h 2 k 2 ] +r 3 



mc . 

(42) 

We introduce the magnetic radius L = y/h/\q\B and de- 
fine £ = k x L + rj q y/L, where r\ q = ±1 is the sign of q. 
Then there is rj q y = £L — k x L 2 and d/dy — (l/L)d/d^. 
The eigenenergies are E n = sE U} k z , where [15| 



E n , kz = v/to 2 c 4 + 2mc 2 huj c (n + 1/2) + (chk z ) 2 . (43) 

The corresponding eigenstates |n) are characterized by 
four quantum numbers: |n) = \n, k x ,k z ,s), where n labels 
the Landau levels, k x and k z are wave vector components 
and s = ±1 label positive and negative energy branches. 
The wave functions are [la] 



* n (r) = (r|n) 



ik x x-\-ik z z / i j "t" 



where </>„(£) are the harmonic oscillator functions 

1 



,(0 



\fZc n 



H„(f)e 



-1/2?" 



(45) 



in which H„ (£) are the Hermite polynomials and C n — 
^2 n nly/jr. We defined fJ^,k x , a = v n,k m ± s/v n ,k x , 
where v n ^ T = ^mc 2 /E n ^ z . 

We want to calculate an average packet velocity in a 
magnetic field. We can, as before, introduce the Heisen- 
berg picture for the time-dependent velocity operator. 
Then the j-th component of the average velocity is, see 
Eq. (HID, 

(fy(t)) = (W\T 3 e lAt/h v J e- im/h \W), (46) 

where Vj = dH/dpj. For the Hamiltonian (|40l) in the 
asymmetric gauge we find 



(r 3 + ir 2 ) 



Px - qBy 



v y = (r 3 + iT 2 ) — , 



(T 3 + IT2) • 

m 



The unity operator is now 

1 = ^ l n )( n K T 3> 



(47) 

(48) 
(49) 

(50) 



where the states (r|n) are given in Eq. (|44| and s n = ±1 
are the quantum numbers associated with the states |n). 
The proof of the above identity is given in Appendix [C] 
Using the unity operator we expand the packet \W) in 
term of the eigenstates of H [see Eq. (Q2 



\W) = J2^\n){n\r 3 \W) = J2^)W n , (51) 



where W n = (n|r 3 |H / ). Inserting \W) into Eq. (PfBl) one 
obtains [see Eq. (pi 



i(t)) =J2 s nSmW*W rn {n\T 3 t 



iHt/h- 



iHt/hi 



(52) 

There is e- iHt / n \n) = e^^n), where w„ = s n E n , k Jh. 
Proceeding the same way as in Section [TT| we have 



{n\r 3 e im ' n 
which finally gives 



(n\e iflU / h T 3 = e^ f (n\r 3 , 



(53) 



(%(*)) = SnS m W r n *W m e<^- w »)*(n|T 3 %|m). (54) 

nm 

The matrix elements of velocity operators calculated be- 
tween the states |n), |m) are 



Xr 



(n|T 3 u y |m 



(nlrs^lm) 



(n|r 3 U2|m) 



-Vn±V m ,k z Sk x ,k'5k z .k' X 



iV2L 

(Vn+ lS m ,n+i - Vn8 m ,n-i), (55) 

V n ,kV rn k z 5 k:E k>8k z ,k' x 



V2L 



(Vn + lS m< n+i + VnS m ,n-i), (56) 

Pz 



'-^n,k z ^m,k z Sk x ,k'Sk z ,k'Sr, 



(57) 



The matrix elements of v y and v x are nonzero for the 
states with m = n ±1 and arbitrary indexes s n and s m . 
The matrix elements of v z are nonzero for m — n and 
arbitrary indexes s n and s m . To simplify the further 
analysis we assume the initial wave packet W(r) to be in 
a separable form 

W(r) = W xy (x,y)W z (z). (58) 

Then there is 

W n = (n\r 3 \W) = ^g z (k z )F n (k z ), (59) 



where 

F n (k x ) = 

in which 



1 



g xy {k x ,y)e-^ 2 R n (0dy, (60) 



Jxy{k x , y) 



/2tt 



w X y{x,y)e lkscX dx, 



(61) 
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and 



9z(k z ) 



For (v y (t)) we obtain 



w z {z)e lk * z dz. (62) 



<%(*)> = c 



E / dfc * 

n «/ — OO 



n,m=0 00 



\g z (k z )\ 2 (vn + l<5 m ,„+i - %Aj<5 mi „_i) [/„, m x 
{ (1 + i^n) cos(w m i - w n i) + (i/^,^ - 1) cos(w m t + cj„t) 
+i{Vm + v 2 ) sin(cj m i - u n t) + i(v 2 n - ^) sin(cj m i + w„i)} . 

In the above expressions we use the notation v n = v n .k z 
and u) n = E nt k z /fi, and 



F*(k x )F m (k x )dk x . 



(63) 



For a Gaussian packet of Eq. (fTTj) one can obtain ana- 
lytical expressions for U„ t m, see Appendix ID] After per- 
forming the summation over m and changing n n + 1 
in <5 mi „_i terms we finally obtain 



A, 



OO 



V n— 

|5*(M| 2 {(^n.+i + ^n) sin(o; n+ it - a; n t) 

- ^n) sin(u;„ +1 i + w„t)} dfc z , (64) 



(«»(*)) 



; ^E 



Vn + 1 ([/„+!,„ + J7„,„+i) 



n=0 

|^(^)| 2 {(1 + ^n+l^n) COs(w n+ ii - U n t) 

+(1 - vl, +1 Vn) cos(w„ + it + w n t)} dk z , (65) 

<«*(*)> = ^E^ / ^i5,fe)i 2 x 

{(1 + i£) + (i _ ^) cos(2w„t)} dfc z . (66) 



Equations (|64|) - (|66|) are our final results for the aver- 
age velocity of wave packet in a magnetic field. Both 
the arguments of sine and cosine functions as well as 
coefficients v n and v n +l depend on k z , so all the inte- 
grals vanish in the limit t — > oo as a consequence of the 
Riemann-Lebegues theorem and the resulting oscillations 
have a transient character. The velocity of the packet os- 
cillates with many frequencies uj n +\ ±w„ (or 2uj n for v z ), 
but in practice the spectrum is limited to a few frequen- 
cies related to the largest coefficients C/ n+ i jT1 and U rhn . 
The frequencies uj n+ i — uj n correspond to the intraband 
transitions and they can be interpreted as the cyclotron 



resonances. These frequencies do not appear in v z veloc- 
ity. On the other hand, the frequencies (d n +i+oj n and 2oj n 
(for v z ) correspond to interband transitions and they can 
be interpreted as the Zittebcwcgung components of the 
motion in analogy to the situation at zero field. The mo- 
tion in the x — y directions requires that k$ x ^ because 
for kjw = all the coefficients U n +i in and U n . n +\ van- 
ish [17J . For the motion in the z direction one needs only 
that k$ z 7^ 0, because the coefficients U n , n are nonzero 
for any ko x vector [l7T |. 

Considering the non-relativistic limit in Eqs. (|B~4l - 
(|66[) . there is hoj c -C mc 2 and hk z -C Trie, so that w„+i — 
uj n ~ hjjj c and ui n +i + u) n ~ 2mc 2 /h. In this limit there 
is Vn+i ~ v n ~ 1, and the ZB part of velocity is nearly 
zero. In this case we may decouple in Eqs. (l6"4")) - (p6")) the 
summation over n and integration over k z . This gives [TtJ 



Vn + lU n+ 

OO 

XX 



n=0 



kp x L 
" V2 



Integrating over k z one gets 



(«»(*)) 

(*x(*)) 

(*«(*)) 



rn 
fokpx 

m 
hk 0z 



sin(w c i), 
cos(ui c t), 



(67) 
(68) 

(69) 
(70) 
(71) 



Thus in the non-relativistic limit the particle moves on 
a circular orbit with the cyclotron frequency in the x — y 
plane and a constant velocity in the z direction. Let us 
introduce a measure of intensity of a magnetic field by 
its relation to an effective Schwinger field KeB s /m = mc 2 
or, equivalently, by L s = K/mc, There is B s = 4.41 x 
10 9 (r7i/m e ) 2 T, where m e is the electron mass. Below we 
perform calculations for pions 7r + having the mass m ~ 
273.1 m e , so the effective Schwinger field is B s — 3.29 x 
10 14 T. 

In Fig. [3] we plot the average packet velocity for three 
values of magnetic field. The ellipsoidal packet is se- 
lected with a nonzero initial momentum ko x . We as- 
sume that the five parameters: d x , d y , d z , L and fc^, 1 
have similar orders of the magnitude which are the op- 
timal conditions for the appearance of Zitterbewegung 
phenomenon. In Fig. [3] we selected parameters: d x = 
0.91(B s /B)X e , d y = 0.82(B S /B)\ C , d z = 0.68(B S /B)X C , 
k 0x — 0.7(B/B S )\~ 1 and k 0z = 0, where B s is the effec- 
tive Schwinger field. For B = 4.5B S we set k$ x = A^ 1 . 
For low fields (B = 0.0045£? s ) the packet moves on a cir- 
cular orbit, see Eqs. (f6"9"|) and ([70]) . For such fields the 
ZB components of the motion are negligible. For higher 
fields the packet motion includes both the intraband and 
interband (ZB) components so that several frequencies 
give significant contributions to the motion. In all cases 
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FIG. 3: Time dependent velocity components for ellipsoidal 
wave packet at various magnetic fields. For low fields, see (a), 
cyclotron motion is obtained, for higher fields, see (b) and (c), 
packet velocity includes both cyclotron and Zitterbewegung 
frequencies. In all cases the motion decays in time. 



FIG. 5: Average packet velocity {v z (t)) in the direction par- 
allel to magnetic field versus time for four values of B. Tran- 
sition to the non-relativistic limit is visible. Parameters are 
the same as those used for Fig. [3] 
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FIG. 4: Average velocity of the spherical wave packet in longer 
time scale. The collapse-revival patterns are seen in ZB oscil- 
lations. The motion has a transient character. 



the motion has a transient character but for low fields 
its decay time is very long. In Fig. @] we plot compo- 
nents of average velocity of a spherical packet in a longer 
time scale. The collapse-and-revival patterns occur for 
both velocity components. After a sufficiently long time 
the oscillations disappear. In Fig. [S] we show the average 



velocity (v z (t)} of an ellipsoidal packet having the same 
parameters as those used in Fig. [3] For large magnetic 
fields the motion in the z direction is similar to that in the 
field-free case exhibiting ZB oscillations, see Fig. [T] For 
smaller fields the ZB oscillations disappear and only the 
classical motion remains, see Eq. ([71) . Finally, it should 
be mentioned that in the two-dimensional case the ZB 
oscillations do not disappear in time [17] . 



V. SIMULATION OF ZB 

The phenomenon of Zitterbewegung for relativistic 
particles in a vacuum has an unfavorable high frequency 
corresponding to the energy gap between the positive and 
negative energy branches: fnoo — 2mc 2 , and a very small 
amplitude on the order of the effective Compton wave- 
length Ar ~ H/(mc), see Eq. (j9|). Thus, similarly to 
the case of relativistic electrons, one can not hope at 
present to observe directly the ZB in a vacuum. However, 
it was recently demonstrated by Gerritsma et al. that 
one can simulate the ZB of electrons in a vacuum using 
trapped ions interacting with laser beams [2| . In this ex- 
periment the authors simulated the linear momentum pi 
appearing in the Dirac equation with the use of Jaynes- 
Cumminngs interaction between the electrons on trapped 
ion levels and the electromagnetic radiation. The deci- 
sive advantage of such a simulation is that one can tailor 
the frequency and amplitude of ZB making them con- 
siderably more favorable than the values for a vacuum. 
Clearly, it would be of interest to simulate the ZB of a 
Klein-Gordon particle using similar methods. The prob- 
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rubber sheet 



string at 
instant t 




•y(x,t) 



element dx 



string's equilibrium position x 



to the wave equation [1 11 ] 



FIG. 6: Classical simulation of KGE according to Morse and 
Feshbach Flexible string is anchored at two points and 
tension T is applied to each end. The string is also attached to 
a thin rubber sheet. At instant t the shape of string is given 
by y(x,t). There are two forces acting on each element dx 
of the string: restoring force Ft due to applied tension and 
elastic force Fk of stretched rubber. 



lem is that in KGE one deals with squares of momentum 
components p 2 , which are more difficult to simulate with 
the Jaynes-Cumminngs interaction. For this reason we 
choose a different route. 

The Klein-Gordon equation appears in several classi- 
cal systems, usually as a modification of the wave equa- 
tion □(/> = 0. Under some conditions KGE is used to 
describe sound waves in ducts 0, HU , electromagnetic 
waves in the ionosphere (2^,[2lJ,_transverse modes of wave 
guides [22| and oceanic waves [23[ . Below we examine in 
more detail a model proposed by Morse and Feshbach 
in which one can simulate KGE with the use of a piano 
string and a thin rubber sheet . Employing this exam- 
ple we demonstrate similarities and differences between 
ZB in the relativistic KGE and its classical analogues. 

Let us consider flexible one dimensional string in the x 
direction, see Fig. [51 We assume that the string is uni- 
form with a linear density p. A uniform tension T is 
applied to each element dx of the string. We neglect all 
other forces acting on the string (e.g. gravity) and the 
stiffness of the string. Let y(x,t) be a displacement of 
the clement dx of the string from its equilibrium position 
at an instant t. We assume that y(x, t) is small com- 
pared to the length of the string and to the distances to 
each end of the string. The restoring force acting on each 
element dx of the string is Ft = Tdx(d 2 y /dx 2 ) and dis- 
placement y(x, t) of the released string changes according 



a 



1 d 2 y 



d 2 y 

dx 2 ' 



(72) 



where u 2 — T j p. Now we attache the string to an elas- 
tic substrate, e.g. to a thin sheet of rubber which can 
shrink or expand in the y direction. Then, in addition 
to the restoring force due to the tension, there will be 
another restoring force due to the elastic rubber acting 
on each element of the string. If the element dx is dis- 
placed to y(x,t) and the rubber sheet obeys the Hook 
law, the restoring force acting on the element dx of the 
string is Fk{x,€) = —Ky(x,t)dx, where K is the elastic 
constant of the rubber sheet. The second Newton law 
for the element dx of the string having mass dm — dxp 
is dxp(d 2 y / dt 2 ) = Ft + Fk, so the equation of motion 
of the released string is 



1 d 2 y d 2 y 
u 2 dt 2 dx 2 



v 2 y, 



(73) 



where v 2 = K/T. Equation ([73")) has the form of wave 
KGE with the light speed replaced by u and the mass 
term m 2 c 2 /h 2 replaced by v 2 . Comparing Eq. ([73]) with 



Eq. ([27]) we find the following correspondence between 
parameters of the two systems 



T 
P 



K m 2 c 



= a: 



(74) 
(75) 



Thus one can simulate values of c and A c by changing 
material parameters p, K and T. 

However, there exist also limitations of such a simu- 
lation and they affect a possibility of observation of ZB 
motion in classical analogues of KGE. The first difference 
between the relativistic KGE and its classical counterpart 
is that the wave function in the relativistic KGE is not 
an observable. On the other hand, all classical analogues 
of </> (such as a displacement of the string, the pressure of 
sound or the oceanic waves, the intensity of electromag- 
netic field etc.) are observable quantities. The second 
difference is that the relativistic function is a function 
of complex variable, while its classical counterpart is a 
function of real variable. A direct consequence of these 
limitations for observation of ZB in classical systems is 
that, for any real function £(r,i) being the solution of 
KGE, the current density associated with this function 
is always zero: j oc [£*V£ — (V£*)£] = 0. Therefore we 
are not able to simulate directly the current or velocity 
oscillations calculated in the previous sections. 

To overcome this problem let us consider the motion 
of a neutral particle described by a real field £. For sim- 
plicity we assume a one-dimensional KGE that can be 
simulated by a flexible string attached to an elastic sub- 
strate described above. In our calculations we use the 
relativistic form of KGE but the final results will be pre- 
sented for parameters corresponding to the flexible string 
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model. We assume the initial wave packet to be a real 
Gaussian function without an initial momentum 



where 



Its Fourier transform is 

w Q (k) = (2d0F) 1/2 cxp[-d 2 fc 2 /2]. 
A real solution £(x, t) of KGE is 
1 f°° 

£(x, t) = — / wo(k) cos(kx — ui k t)dk, 
2n .1^ 



(76) 



(77) 



(78) 



where ui k = ujq y/l + (k\ c ) 2 . The average current for the 
real wave packet of Eq. (|76p is zero and no ZB occurs. 
Thus we turn to other physical operators which do not 
commute with the KG Hamiltonian (JlJ . Namely, we cal- 
culate a variance of the position operator for the above 
real function £(x, t) 



(79) 



since {<fi\x\(f>) = 0. Assuming £(x,t) in the form (f78|) we 
have 



V = Jjj wo(k)wo(k') cos(kx — ui k t) x 
cos(fc'x — Lu k d)x 2 dxdkdk' = 

B k B k ,e lx{k+k) + B k B* k ,e lx{k - k,) 

BlB k ,e- tx{k - k ' ] + B* k Bl,e- ix{ - k+k 'A x 2 dxdkdk' ,(80) 



where B k = wo(k)e ~ lUkt / (4ir). Consider the first of 
the four terms given above. Because Wo(k) and B k de- 
cay exponentially for k — ¥ ±oo, one can change the or- 
der of integration over x, k and k' and replace x 2 — > 
(d / dik)(d / dik'). Then we integrate by parts over k 
and k' and obtain 



B k B k ,e lx{k+k) x 2 dxdkdk' 
dB k 9B k ' elx(k+k)dxdkdkl 

dk. 



= 2tt 



dk dk 1 

dB k 8B k , 



dk dk' 



(81) 



The other three terms in Eq. ([5U|) are calculated similarly. 
After some manipulations we find 



V = Vf + v° sc + V 2 C + v° sc + v 3 , 



(82) 



Vf = 



e- dk '(kdfdk, 



V° 



e~ a k (kd) z cos(2uj k t)dk, (84) 



vs = 



-rrOSC 
V 2 



Vo = 



d 3 

2^ 
d 3 

2^ 

djctf 

^ Joo 1 + (£A C ) 2 

\kX c 



(83) 



2 /-oo „-d z k' 



{k\ c f 



dk, 



(85) 



2 /-oo „-d'k z 



d(ctf 

7oo l + (k\c) 2 

d 2 (ct) f°° e- d2k \k\ c ) 
V 1 + ( kX c) 2 



■cos(2cj fc t)dfc,(86) 



sin(2uj k t)dk. (87) 



The term V3 is odd in k and it vanishes upon the in- 
tegration. For t — the variance in Eq. (|82|) is equal 
to the variance Vb = d 2 /2 of the initial packet wo(x). 
The variance given in Eq. (|82l) consists of oscillating and 
non-oscillating terms. For large times the non-oscillating 
terms grow in time as d 2 /2 + Ct 2 , where C is a con- 
stant depending on d. The quadratic dependence of the 
variance on time is similar to that of a Gaussian wave 
packet in the non-relativistic quantum mechanics. The 
oscillations in Eqs. (l84l and ((86)) have the same interband 
frequency 2u> k as the velocity oscillation given in Eq. . 
Therefore the oscillations of variance of the position oper- 
ator can be interpreted as a signature of Zitterbewegung 
in classical systems. The term V° sc has a decaying char- 
acter and it vanishes after a few oscillations. The V£ sc 
term gives persistent oscillations because of the presence 
of t 2 factor in front of the integral. To estimate the time 
dependence of these oscillations we consider the limit of 
large packet widths d ^S> X c . In this case the Gaussian 
function in Eq. (|56"j) restricts the integration to small val- 
ues of k. Then we may disregard (fcA c ) 2 term in the de- 
nominator of integrand and expand u k under the cosine 
function. This gives approximately 



V° sc 



d(ct) 2 
' 2^ 
d{ctf 



E 

r)=± 



l (d 2 + ii]uj t) 3 / 2 



For large time we may approximate in Eq. 



(kX c ) 2 cos[o;oi(2 + k 2 X 2 )]dk 

2irjLdot 

(88) 



-C d t 1/2 cos(2w i), 



(89) 



where Cd is a constant depending on d. Thus the os- 
cillations of variance are persistent, their amplitude in- 
creases with time as t 1 / 2 and their frequency is 2ujq. Since 
non-oscillating terms V£ increase as t 2 , the total variance 
of the packet has a quadratic time dependence with su- 
perimposed oscillations. This behavior is illustrated in 
Fig. [71 In our classical considerations we do not face the 
problem of negative variances that can occur for some 
quantum systems, see Refs. [U [HI- For t < 5i c the 
oscillations have an irregular character because of the 
contribution of V° sc term. 
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FIG. 7: Calculated classical variance of position of wave 
packet propagating according to KGE. Dashed line: non- 
oscillating part of variance; solid line: total variance. For 
string oscillations analyzed in the text there is = 4.47 mm 



and ti 



2.37 x 10" 



Estimating the characteristic frequency 2wo for the 
flexible string attached to elastic substrate we have 



w 2 



1 A 

m c 



c 2 x 



A 2 . 



K 
P 



(90) 



so that the analogue of the relativistic frequency ujq does 
not depend on the applied tension. Taking a piano cop- 
per string of the bulk density p^D = 8940 kg/m 3 and 
having cross section of radius r = 1 mm one gets a linear 
density p = nr 2 p 3D = 2.81 x 10~ 2 kg/m. We identify 
the rubber elastic constant K with the Young modu- 
lus K = 0.05 x 10 9 N/m 2 . Then the analogue of ZB 



frequency given in Eq. ([9"0]) is 2wq 



.44 x 10 4 s"\ i.e. 



the corresponding frequency is fo = 13.43 kHz, which 
can be heard by the human ear. The characteristic time 
of ZB oscillations is t s c = l/wg = 2.37 x 10~ 5 s. As- 
suming the tension of the string T = 1000 N we find 
from Eq. (175)) that the simulated Compton wavelength 
is A* = 4.47 mm. The initial wave packet should have 
widths d on the order of a few A*, i.e. of a few centime- 
ters, and it will move with the velocity u = 188.7 m/s, 
see Eq. (f74"|) . Thus, it is really possible to simulate and 
observe the Zitterbewegung phenomenon in this system. 
Finally, we observe that in classical simulations all the 
involved quantities are well defined observables. Since 
classical KGE does not reproduce but only simulates 
the quantum KGE, we are allowed to consider quanti- 
ties which are not well defined in the quantum world. 



VI. DISCUSSION 

Our main results for ZB of KG particles in absence 
of fields are shown in Fig. Q] and in the presence of a 
magnetic field in Figs. [3] - [3] It is not our purpose here 
to consider difficulties of the one-particle Klein-Gordon 
equation but we keep them in mind. In particular, we 
do not consider particle tra ject ories as they are believed 
to be not well defined, see On the other hand, we 
describe average particle velocities and currents both in 
the Hamiltonian and wave formalisms. The results can be 
compared to those for relativistic electrons in a vacuum 
described by the Dirac equation as well as for electrons 
in solids. 

Similarly to the Dirac electrons, the ZB phenomenon 
of KG particles is due to the interference of positive and 
negative energy states. In the non-relativistic limit one 
of the two components progressively vanishes and the 
ZB contribution to the motion disappears. This can be 
clearly seen in Figs. [3] and [5] as well as in Fig. 3 of Ref. [TtJ 
for the Dirac electrons. If particles are described by wave 
packets the ZB motion decays in time, see our Fig. [T] for 
KE particles and Fig. 2 of Ref. 26] for the Dirac elec- 
trons. This is a general consequence of the Riemann- 
Lesbegues theorem, as indicated by Lock [3] , calculated 
by the present authors [27j and experimentally confirmed 
by Gerritsma et al. • In all cases the basic frequency of 
ZB oscillations is given by the energy difference between 
the positive and negative energy branches: Tuujz — 2mc 2 
with the corresponding particle mass. The main differ- 
ence with the Dirac electrons is the spin. For KG parti- 
cles the interband ZB frequencies in a magnetic field do 
not include the spin energies, one does not deal with the 
Fermi sea for the negative energy branches, etc. The 
KG Hamiltonian is quadratic in momenta which does 
not allow a direct simulation with the use of Jaynes- 
Cumminngs interaction. 

As to the Zitterbewegung of electrons in narrow-gap 
semiconductors and in particular in zero-gap monolayer 
graphene, one should emphasize that, although it is also 
described using a two-band model of band-structure [3] , 
its physical nature is completely different from ZB of 
particles in a vacuum. The ZB in semiconductors or 
in graphene results from the electron motion in a peri- 
odic potential [28j . In zero-gap situation in graphene the 
ZB frequency is given by the difference of energies be- 
tween positive and negative energy bands corresponding 
to the average value of quasi- momentum hko for the wave 
packet [27]. A one-dimensional system which strongly 
resembles the KG particle in a vacuum is presented by 
electrons in carbon nanotubes: one can neglect the elec- 
tron spin dealing with an energy gap controlled by the 
tube's diameter (29|. The resulting ZB frequency and 
amplitude have values easily accessible experimentally. 
On the other hand, it is at present not clear how to fol- 
low dynamics of a single electron in a solid. As to KG 
particles in a vacuum, one is bound to recourse to simu- 
lations since the ZB frequency and amplitude as well as 
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field intensities necessary to see ZB effects in the pres- 
ence of a magnetic field, exceed the present experimental 
possibilities. 

We present a classical simulation of ZB by using a me- 
chanical system and calculate the oscillating variance of 
position of the wave packet. The variance of position op- 
erator for the Dirac Hamiltonian was calculated by Barut 
and Malin 30] who found it to be the reminiscence of ZB 
of electrons in a vacuum. The present authors analyzed 
in Ref. [13] the variance of position operator in bilayer 
graphene and found its oscillating character with the fre- 
quency equal to that of ZB. 

One should finally remark that the attempts are con- 
stantly made in the literature to overcome the above men- 
tioned difficulties in the interpretation of position oper- 
ator in KG equation. In particular, Mostafazadeh [3l[ 
proposed a redefinition of the scalar product of solutions 
to KGE which allows one to obtain positively defined 
probability distribution of position. Semenov et al. [25| 
proposed to limit the allowed solutions of KG equation to 
those having positive-definite probability distributions. 
They showed that the physical solutions of KGE fulfill 
this criterion. If the above attempts are accepted one 
could analyze ZB of the position operator for KG parti- 
cles, see Eq. ©. 



VII. SUMMARY 

We considered the trembling motion (Zitterbewegung) 
of relativistic spin-zero particles in absence of fields and 
in the presence of a magnetic field using the Klein- 
Gordon equation. We aimed to describe physical observ- 
ables (currents and velocities) calculating quantities aver- 
aged with the use of Gaussian wave packets. Surprisingly, 
the calculated particle velocities can exceed the veloc- 
ity of light for sufficiently large momenta indicating that 
KGE docs not posses an automatic restriction of relativ- 
ity. We showed that the trembling motion has a decaying 
character resulting from an interference of positive and 
negative energy sub-packets moving with different veloc- 
ities. In the presence of a magnetic field there exist many 
intcrband frequencies that contribute to Zitterbewegung. 
On the other hand, in the limit of non-relativistic energies 
the interband ZB components vanish while the intraband 
components reduce to the cyclotron motion with a single 
frequency. The trembling motion was simulated using 
the classical system obeying the Klein- Gordon equation 
- a stretched string attached to a rubber sheet. The cal- 
culated variance of position of the sting shaped initially 
as a Gaussian packet exhibits oscillations corresponding 
to Zitterbewegung with the correct frequency. 



field. We define the creation and annihilation operators 

(Al) 



a = (S + d/d£)/V2, 
a+ = (£-0/00/ A 



and rewrite Eq. (|40[) in the form 



+ T 3 mc 



(A2) 



where w c = qB/m is the cyclotron frequency and T 
(T3 + vr%). The current density is 



2im 

e 



A^T 3 t ip, 



(A3) 



and the average current is (j) = j jd 3 r. We introduce 
the current operator J in such way that for j given in 
Eq. (|A3]l there is 



Ml 



jd 3 r. 



(A4) 



Note the presence of T3 in the matrix element. In the 
asymmetric gauge one has 



mc 



(A5) 



2m 



dip dip^ 



(h) = -7T- / l^T^ - -£-T 3 Ti>) d*r.(A6) 



dy dy 



Below we assume the function if> to be Gaussian-like. In 
that case we may simplify the above expressions for the 
average current by integrating by parts the terms includ- 
ing derivatives of ip' 



Ox) 

(jy) 



ih 

m 

ih 

rn 



■I t- I <l 'v — 

ox J mc 



9 ^ ^-^lU^fyAm 
mc J \ / 

/ {^i) " 3r ' (A8) 



so the components of the current operator are 

- ih a- d qB , . „. 

J* = T^-—Ty, (A9) 



J y ^ -- f 7T- ( A1 °) 

m dy 

In the Heisenberg picture the time-dependent current op- 
erator is 



Appendix A 



J{t) = e i6t ' h J{Q)e- iAt ' 1 



(All) 



In this Appendix we calculate an exact time depen- 
dence of current operators for a KG particle in a magnetic 



where H is given in Eq. (|A2[) . Our task is to calculate the 
time evolution of the current operator J x {t) and J y (t). 
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By averaging these operators over the state ip, as shown 
in Eq. (|A4[) , one obtains the time-dependent charge cur- 
rent corresponding to ip. 

It is convenient rewrite current operators in Eqs. (|A9|) 
and (|A10p in the form 



J, = V- 



qB 



m y/2mc 
ih 



m 



where we introduce three auxiliary operators: 

P = ^ + i ^lx^ 
J = (r 3 + iT 2 )a = Ta, 
J + = (T 3 +iT 2 )a+ = fa+. 



(j + J + ) , (A12) 
J-J + ), (A13) 



(A14) 

(A15) 
(A16) 



We calculate the time dependence of J , j + and V in 
a way similar to that described in Ref. [17]. Consider 
first the operator V . From the equation of motion Vt — 
(i/h)[H,V] one has 

V t = — [n,V] = 2twon^, (A17) 

where we used T 2 = 0. Since {H, Vt} = 0, there 
is [H,V t ] = 2HV t - {H,V t } = 2HV U and one obtains 



Vtt = %HV t . 



(A18) 



We solve this equation for Vt and then integrate the so- 
lution over time 



V(t) = —re 2iHt/h Vt(0)+C, 
2iH 



(A19) 



where C is a constant of integration. Applying the initial 
conditions: fi(0) = f(d/dx), V t {0) = 2iu; n(d / 'dx) , and 
using the identity H^ 1 = H/E 2 we have 



( e 2Wt/h _ ^ 



d 



l)n-. (A20) 



It is seen that V(t) in Eq. (|A20[) satisfies the initial con- 
ditions for V(0) and Vt(0). The form oiV(t) given above 
resembles results obtained for the position operator in 
the field- free case by Fuda and Furlani 

Now we turn to the operators J and J + . From 
Eqs. (|A15j) and (|A16|) one has 



J t + = 2zw Tia + , 
where loq = mc 2 /h. Using [a,a + ] = 1 one obtains 



{H,J t } 
{Hj t + } 



-2iuJohijj c j + . 



(A21) 
(A22) 

(A23) 
(A24) 



Upon applying the identities 

[H,J t ] = +2HJ t -{H,J t }, 
[Hj t + ] = -2J+H+{H,J+}, 



we get 



J tt = +(2i/h)HJ t -2uj u c J, 
J+ = -(2i/H)J t + H-2u> uJ + 



(A25) 
(A26) 



(A27) 
(A28) 



In Eqs. (|A27I) and (|A28|) we eliminate terms with the first 
derivatives using the substitutions J = exp(+iHt/H)B 
and J + = B + exp(— iHt/H), respectively. This gives 

Btt = -(n 2 + 2uj c u;o)B, (A29) 
B+ = -B + ((l 2 +2cj c w ), (A30) 

where Cl — H/h. In the above equations the operator 



M 2 = fi 2 + 2w c w 



(A31) 



stands on the left-hand side of B, but on the right-hand 
side of B + . Solutions to Eqs. (|A"29)) and (|A"30]) are 



B 
B+ 



iMt f 



JMt/ 



+ AMt 



(A32) 
(A33) 



where M. = +V M 2 is the positive root of jCi 2 . 
Both Ci and C 2 and their complex conjugates are time- 
independent operators. 

Using the initial conditions: B(0) = j(0) — Ta and 
B t (0) = J t (0) = +2iu> Tia [see Eq. (|A21|) ] and similar 
expressions for B + (0) and Bf (0), we find that J{t) = 
Ji(t)+Mt), where 



J x (t) = l e Mt e -iMt 



± JQ.t +iMt 

2 



i(0) + M- l J{Q)h ,(A34) 
i(0) - M- l J(ti)ti\ .(A35) 



Similarly, one can express J + {t) = (t) + J 2 + (t) , where 



e +iMt e- iU tA36) 



e -iMt e -in 



(A37) 



The results are given in terms of the operators f2 and M. . 
To finalize the description, one needs to specify the phys- 
ical sense of functions appearing in Eqs. (IA34[) - (IA37|) . 

For a reasonable function f(D) of an operator Z), 
its eigenenergies Ad and its eigenstates |d), there exists 
the following relationship: /(Z?)|d) = /(Ad)|d), provided 
that /(Ad) exists. To find the meanings of the opera- 
tors M" 1 and e ±lMT we express them as functions of 
the operator M 2 = H 2 /h 2 + 2uj uj c , see Eq. (|A3L|) . From 
the definition of M 2 it follows that its eigenstates are 
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equal to the eigenstates |n) of H . The eigenvalues A 2 of 
the operator M 2 are A 2 



p2 



and we obtain 



M^n) = (M 2 ) ±1 / 2 \n)=r 1 E^ k _\n), (A38) 



e ±iMt \n) = e ±i ^ M y /t \n)=e Mr ' E -+ 



where 77 = +1 or r\ = —1. As seen from Eqs. (|A34|) - 
(|Al7)l . the sums Ji(t)+J 2 (t) and J+ (t) + J+ (t) do not 
depend on the sign of 77, so we select r\ = +1. 

Finally we show that the matrix elements of the op- 
erator j{t) = j\{t) + Ji{£) are equal to the matrix el- 
ements of the current operator Jn(t) — e int J(0)e~ lflt 
in the Heisenberg picture. The operator J is propor- 
tional to the annihilation operator a whose non- vanishing 
matrix elements are (n'\a\n) = \/n~+T5 n ' , n +i, so we 
select two eigenstates of KG Hamiltonian |n) = \n,s) 
and |n') = \n + l,z), see Eq. (l4"4"l) . Here we omitted 
quantum numbers k x and k z . For Jnii) one has 



(u\T 3 J H (t)\n') 



glSUJnt g — lZUJn+lt 



lt J{0)r 



consider (1,1) component of the velocity operator for a 
KG particle given in Eq. ([7]). For the wave packet (r\w) = 
w(r)(l, 0) T with one nonzero component the average ve- 
locity is given by the average of (v z )u(t) over the func- 
tion w(r). The unexpected feature of operator (v z )n(t) 
is that for large p this velocity can exceed the speed of 
light c. 

There are two possible ways to overcome this problem. 
We can additionally assume that \p\ < mc, which ensures 
that the velocity (v z )n(t) does not exceed c. This con- 
dition is equivalent to \q\ < 1 in the text, see Eq. (J7]). 
Alternatively, one can take the initial wave packet w(r) 
which does not contain components with \p\ > mc. Then 
the Gaussian packet in Eq. (fT2|) must be replaced by a 
non-Gaussian packet w'(r) of the form 

(k\w') = (2dy^) 3/2 exp[-rf 2 (fc - fc ) 2 /2]e(A c - |fc|), 

(Bl) 

where 0(£) is the step function. 

For the Dirac Hamiltonian Hd = cj\. ttjPj + mc ft, 
(A40) the situation is different. Expanding e lH ot/n m a p 0we r 



where we define l 7(0) nn ' = (n|T 3l 7(0)|n'). To calculate 
the matrix elements of j\[t) we use Eqs. (IA38|) - (|A39[) 
and obtain 



(n\M- l J{Q)Q,\v!) 



which finally gives 



E 



n+l 



zE n+1 



Zj(0)r 



(A41) 



(n\Mt)W) 



l + z 



l 7(0)„n'e IS "" t e-^"+ lt , (A42) 



1 - z 



(n|J 2 (t)|n') = ^ 1 7(0) na /e <au »*e +fa '»+ l *. (A43) 

The matrix elements of J\ (t) are nonzero for z — +1 only, 
while the matrix elements of J% (t) are nonzero for z = — 1 
only. Comparing Eqs. (|A42|) and (IA43|) with Eq. (IA40j) 
we see that for each of four combinations of s = ±1 
and z = ±1 the matrix elements of J7jy(i) are equal to 
the matrix elements of j{t) = j\{t) + jiif), which is 
what we wanted to show. Calculations for J + {t) are 
similar to those for J{t). The compact equations (|A34j) 
- (| A37|) are our final results for the time dependence 
of j(t) and J + {t) operators. These equations are ex- 
act and they are quite fundamental for relativistic spin- 
particles in a magnetic field. If we calculate average 
currents of Eqs. (|A12[) and (|A13[) with the use of ex- 
pressions (IA42P - (|A43|) and the wave packet (fTTj) . one 
obtains results corresponding to the velocities given in 
Section EH 



Appendix B 

In this Appendix we analyze in more detail the rela- 
tion of the particle velocity to the speed of light. We 



series one obtains an expression analogous to e 



iHt/h 



in Eq. 



After some algebra we find 



(««)£(*) 



mc 2 p z 
m 2 c 2 + p 2 



[1 - cos(2Et/H)] 



given 



(B2) 



In contrast to the KG case, the velocity operator given 
in Eq. (|B2|) has correct relativistic behavior for all val- 
ues of p. In Eq. (|B2I) the expression in square brack- 
ets oscillates between zero and two. The factor v D (p) = 
mc 2 p z / (m 2 c 2 +p 2 ) tends to zero for both large and small 



values of p. Its maximum is at p 7 ' 
which one obtains 



(». 



>?,(*) = ! 



1 



cos 



(2V21 



(0,0, mc) for 



(B3) 



The above velocity never exceeds the speed of light. 
Therefore, when calculating the average velocity of the 
wave packet for the Dirac Hamiltonian, there is no need 
for an artificial truncation of the high momentum com- 
ponents of the wave packet, as proposed in Eq. (|B1[) for 
a KG particle. 



Appendix C 

We prove here some identities appearing in the pre- 
vious sections. We begin with the identity in Eq. (|50[), 
Closing Eq. ([50|) with the use of states (r\ and \r'), em- 
ploying Eq. (|44|) and writing explicitly the summations 
and integrations over the quantum numbers we obtain 



n 



n=0 



■DC 



x 1 e ik * {x ~ x ' ] dk x I e^-^dkz 

-OO J —OO 



IX. 1 



E 

s=±l 



ST 3 , 



(CI) 
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where /i* = ^ k s In the above equation the sum- 
mation over n gives 5te>, the product of the two inte- 
grals is 4:7T 2 S XtX iS z ^ Z ' , so the product of the three terms 
equals 4:Tr 2 S r ^ r i. Taking the explicit form of ^ =v±sjv 
where v — -J 'mc? j 'E n we obtain for the last line of 
Eq. (JCU) 



/ s(v + s/ v ) 2 s(v 2 -l/v 2 ) 
2^ I s{v 2 - \jv 2 ) s{v-sjv) 2 

s=±l 



1 

-1 



(C2) 

Collecting all numerical factors we see that the right hand 
side of Eq. (IC1|) equals b~ r y- 

Next we prove the identity used in the derivation of 
Eqs. (fl"9|) and ([53]) . Let O be any operator for which O = 
t 3 &t 3 , where dagger signifies the Hermitian conjugate. 
We want to show that 



Eqs. (I64|) - (|66|) . Here we assume the initial wave vec- 
tor in the form fe = (ko x , 0, ko z ). Using the definitions 
of g xy (k x ,y), F n (k x ) and U m ,m we obtain (see Ref. H3) 



9xy (k x ,y) 



ird. 



(Dl) 



and 



F n (k x ) = '^ U * e -h<(^-^)\~^ Un (~k x c), 



^/2ird y C n 



(D2) 



where D = L 2 / J L 2 + d 2 y) c = L 3 /JL 4 - d y , and 



T 3 e a = e at T 3 . (C3) 
To this end we expand the exponents and get 



T3 



1 



6 o 2 

IT + ~2l 



1 



0t ota 



1! 



2! 



r 3 . (C4) 



Since = r 3 O t T 3 , there is £>t = T 3 e>T 3 . Then for n > 
there is O tn = r 3 0™T 3 and we obtain for the RHS of 
Eq. El 



Ot 2 



1 ot 

1 + T! ■ 2! 



7-3 



e> 2 

1 + T 3yy T 3 + r 3-^r r 3 



T 3 



r 3 



O 

IT 



e> 2 

~2j 



r 3 e 



o 



which is the desired result. 



2 j2 \ " /2 



-4, 



!2ivdy I L — d y 



L 2 + d 2 \ L2 + d l 



(D3) 



TT _ A m A n LQd x ^n e \ - /../ /;/ \/ // 



Z=0 



x((l - (cQ) J H m+ „_ 2 / \ -j==== I , (D4) 



in which Q = l/^d 2 + D 2 , W = d x DQk 0x , and Y = 
d 2 ko x Q. For the special case of d y = L, the formula 
for U mn is much simpler: 



exp 



C m C n L 
1 2 k 2 . 



2P 



m+n+1 



2P 2 



H 



rn+n 



. (D5) 



Appendix D 

In this appendix we quote for completeness all for- 
mulas necessary for a calculation of coefficients U mn in 



where P = y d 2 + \L 2 ■ In the above expressions the 
coefficients U mn are real numbers and they are symmet- 
ric in m, n indices. For further discussion of of U m ,n see 
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